Exponential growth of eccentricity in secular theory 
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The Kozai mechanism for exponentially exciting eccentricity of a Keplerian orbit by a distant 
perturber is extended to a general perturbing potential. In particular, the case of an axisymmetric 
potential is solved analytically. The analysis is applied to orbits around an oblate central object 
with a distant perturber. If the equatorial plane of the central object is aligned with the orbit of 
the distant perturber (axisymmetric potential), a single instability zone, in which eccentricity grows 
exponentially, is found between two critical inclinations; if misaligned (non-axisymmetric potential), 
a rich set of critical inclinations separating stable and unstable zones is obtained [J]. The analysis 
is also applied to a general quadratic potential. Similarly, for non-axisymmetric cases, multiple 
stability and instability zones are obtained. Here eccentricity can reach very high values in the 
instability zones even when the potential's deviation from axisymmetry is small. 



A striking aspect of the secular evolution of a Keple- 
rian orbit weakly perturbed by a distant orbiting mass 
is that the eccentricity can grow exponentially to high 
values @, Q . This so-called Kozai mechanism operates in 
a finite range of mutual inclinations 39.2° < i < 140.8°. 
The high eccentricity excited from initially nearly circu- 
lar orbits by the Kozai mechanism is suggested to play an 
important role in the formation and evolution of many 
astrophysical systems [e.g. - In this letter we general- 
ize this mechanism by studying the stability of precessing 
circular orbits for a general perturbing potential. 

Consider a test particle orbiting a central mass M sub- 
ject to a small, time-independent perturbation by a po- 
tential — $(r), so that the total potential is 



V(r) 



GM 



*(r) 



(1) 



At any given time the orbit is approximately Keplerian, 
with the orbital parameters changing with time. In the 
secular approximation, the Hamiltonian is averaged over 
the orbit's phase (mean anomaly) to obtain equations of 
motion for the orbital parameters. Under this approxi- 
mation, the semi-major axis a is constant in time. The 
problem reduces to understanding the long-term evolu- 
tion of the remaining 4 orbital elements describing an 
orbit at the given semi-major axis. 

It it useful to describe the orbi t of th e particle by two 
dimensionless vectors: j = J/yGMa, where J is the 
specific angular momentum vector; e, a vector pointing in 
the direction of the pericenter (point of closest approach) 
with a magnitude equal to the eccentricity e. Note that 
e and j satisfy 



.1 



|j| 2 = l- e 2 , 



e-j =0, 



(2) 



leaving 4 independent parameters. 

The equations of motions for these variables are 



(It 



= j x Vj0 + e x V c 
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= jxV e ^ + exV j( )i, (3) 
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where 



ffl,e) = ^-, r = t/t 
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and ($) (j, e) is the potential time-averaged over an orbit 
set by e, j and the fixed a, which is conserved in time. 
When using equations ([3} , the 6 components of e and j 
should be considered independent variables, while only 
the solutions that satisfy the physical conditions Eq. @ 
should be considered. Note that these physical condi- 
tions represent a gauge freedom that does not affect the 
equations of motion [9j . 

Throughout this letter the following examples are con- 
sidered. 

1. "Kozai" - The perturbing potential is produced by 
a distant orbiting mass M per at semimajor axis a por 3> 
a. The leading term in the expansion in powers of 
a/ctper (quadrupole=(a/a pcr ) 2 ) of the potential averaged 
in time over the perturber's orbit, is given by $Koz(r) = 
GMp Cr (463 cr )' 1 [r 2 -32; 2 ] where 6 por = a pcr (l-e 2 CI .) 1/2 is 
the semi-minor axis of the perturber's orbit and z is the 
direction of the perturber angular momentum. By choos- 
3GM pcr a 2 /(86 pcr 3 ), the resulting potential 



ing $ ,Koz 

in terms of e and j is Q 



'Koz 



Jz 



5e 2 + 2e 2 



(5) 



2. "Oblateness" - The perturbing potential is the 
quadrupole potential arised from the oblate central body 
and is given by $ bi(r) = GMJ 2 R 2 (2r 5 )- 1 [r 2 - 3z 2 ], 
where 3i is the gravitational quadrupole coefficient and z 
is the direction of the central body's spin axis. By choos- 
ing $o obi = 3J2GMR 2 / (16a 3 ), the resulting potential is 
1 



^Obl 



= 4 



Jz 



3 J 



r 



(6) 



3. "Koz-Obl" - The perturbing potential is a combi- 
nation of the above two and is given by 



'Koz-Obl 



(7) 
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where $o,Koz-Obi = ^o.Koz is chosen and 
e = ( 5o,Obi/ < £'o,Koz- Note that the z axes of <j>Koz and 
4>Obi are not necessarily the same and the inclination 
between them is denoted as «koz-om- 

4. "Quadratic" - The perturbing potential is a general 
quadratic function of the spatial coordinates, and by a 
proper choice of coordinate system, it can be expressed 

as$ Quadr atic(r) = (2 7 rG/3) / 9 cff [(l + ( 5 + ( T)r 2 -3(z 2 + V)]- 
The parameter p e g has dimension of density, and is re- 
lated to the local matter density pi oca \ by the Poisson 
equation, which reads p e ga = pi oca i- Such a quadratic 
potential represents the leading order of the potential 
that arises from any mass distribution that is approxi- 
mately constant in the vicinity of the test particle's or- 
bit. In particular, it reduces to the Kozai potential when 
8 = a = 0. Examples include the Galactic tide [e.g. 
4] and a combination of multiple distant orbiting per- 
turbers. By choosing $ , Quadratic = nGp c sa 2 , we get 

Quadratic = j 2 z + 5j 2 - 5(e 2 z + 6e 2 y ) + 2(1 + 5 + a /2)e 2 . (8) 

Formulation of the Problem The problem under study is 
to find the conditions under which Eqs. Q have process- 
ing circular orbit solutions that are unstable and lead to 
exponential growth of the eccentricity. 

For circular orbits to be a solution to Eqs. (|3J|, it 
is sufficient that V e <^|e=o = 0. Henceforth we consider 
potentials with reflection symmetry $(r) = <£(— r) for 
which the leading order in e is quadratic and thus they 
satisfy the above requirement. To the first order in e, 
equations Eq. © can then be written as, 



dr 
de 



j x Vj<? 
■ M(j)e, 



(9) 
(10) 



where the coefficients of the matrix M(j) de- 
fined according to de^ / 'dr — M^ n e n are given by, 

Mkn = £klm3lde m ,e n <P\e=0 + £knmdj m 0| e =O ■ 

Circular orbits precess periodically according to Eq. 
©. To see that the precession is periodic, note that j 
is confined to the sphere |j| = 1 and follows the closed 
contour lines of </>(j,e = 0) =const. 

To first order in e, the trajectory of j(t) is the same as 
that of the circular orbit. The problem reduces to solving 
the linear equation Eq. (|10p , with time- varying, periodic 
coefficients according to the periodic trajectory of j(t). 

Axisymmetric potentials Consider first axisymmetric 
potentials (with reflection symmetry). These include the 
Kozai, oblateness, quadratic (with (5 = 0) potentials in- 
troduced above as well as other potentials such as the 
averaged potential of a perturber on a circular orbit to 
all orders in the multipole expansion Q or the potential 
of a proto-planetary disk [l2[ • By eliminating j using the 
gauge freedom j 2 = 1 — e 2 , such potentials can be written 
in the following form 



where ao, a/ and a z are functions of j z . Note that 
V e </> = aje + a z e z z. In this case the unperturbed or- 
bit's angular momentum precesses around z according 
to dj/dr — tlz x j (Eq.©), where tl = dfl/dr = —dj^<j) 
is the constant precession frequency and Q is the lon- 
gitude of the ascending node (angle between z x j and 
x). It is useful to work in a rotating reference frame, 
z'=j, x = z x j / sini and y = z x x processing with 
j, where i is the orbit's inclination (angle between j and 
z) which is constant in the linear approximation. Note 
that e x i = e cos uj,e y i = esinw where lo is the argument 
of pericenter (angle between e and z x j). 

The linear equation, Eq. (llO[) . written in this frame, and 
restricted to physical eccentricity vectors perpendicular 
to j (vectors satisfying e z i = 0) is 



d 



e y i 



—a/ — a z sin i 
a! 



(12) 



In Eq. ()T2|) , the coefficients a/ and a z are to be evaluated 
at e = and thus are functions of i alone and are time 
independent. The eigenvalues are A± = ±\/A, where 



-ai(ai 



(13) 



and eccentricity grows exponentially if and only if A > 0. 
For the cases of Kozai, oblateness and quadratic with 
6 = (Eqs. ©, © and ©), a/, K oz = 4 and 
a z ,Koz = -10, a/.obl = 20 cos 2 i - 4 and a Zj0 bi = 0, 
and a/quadratic = 4 + 2cr and a Z: Q uadr atic = -10, respec- 
tively. The instability criterion for the Kozai Mechanism, 
i > i c = sin _1 (-\/2/5) @, Q, is reproduced while circu- 
lar orbits are stable for all inclinations in the oblateness 
case. 

Consider next a test particle perturbed by the com- 
bined effects of oblateness and Kozai (Eq. (0) in the 
case where they are aligned (i.e. iKoz-ow = 0). One ex- 
ample is a satellite of Jupiter perturbed by its oblateness 
and the Sun (the spin axis of Jupiter is aligned with its 
orbital angular momentum to about 3°). In this case (f> 
has coefficients a] = 4[1 + e(5cos 2 i — 1)] and a z = —10. 

Using the instability condition A > 0, for e < l(e > 1) 
instability occurs at inclinations i > i c \ (in the range 
i c .i < i < i c< 2) where [l| 



sin i c i 
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1 



4e + 2 
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(14) 



4> = a o + T;aie 2 + -a z e 2 + 0(e 4 ), 



(11) 



In the limit e — > (oo), the instability zone of Kozai 
(oblateness) is reproduced. Note that for large e the 
unstable zone reduces to a small interval with a width 
Ai « 2/ (5e) in the vicinity of the so-called "critical in- 
clination" i c = sin _1 ((4/5) 1 / 2 ) [see e.g. HH for discussion 
of the critical inclination]. 

General perturbing potentials Consider next a general 
perturbing potential. The analysis can be proceeded sim- 
ilarly to the axisymmetric case by working in the above- 
mentioned rotating frame, where in general i is not con- 
stant but varies periodically with time. The equations of 



3 



motion (flOl) have the form de' jdr = m(t)e' where m is 
a 2 x 2 matrix with coefficients which change with time 
(periodically) and therefore cannot be solved analytically 
in general. 

The stability of the solution e' = (precessing circu- 
lar orbit) can be studied by studying the linear operator 
u which represents the integration of the above equa- 
tion over one period T and is defined by e'(T) = ue'(O) 
[Floquet anlysis, e.g. [H, section 28]. After n periods 
(t = nT), e'(nT) = u n e'(0). The condition for instabil- 
ity is that u has an eigenvalue x with magnitude larger 
than unity. The corresponding eigenvector e/ satisfies 
u"e/ = x n ej. Given that \x\ > 1, any initial condition 
results in exponential growth of e with time (except pos- 
sibly for initial e which is perpendicular to ej). 

Given that e'(0) = e(0) and e'(T) = e(T) (since the 
rotating frame returns to itself after one period), the 
linear operator u can be calculated using Eq. (fTUj) di- 
rectly, without moving to the rotating frame. This can 
be done by choosing any two physical initial conditions 
e;(0), (I = 1, 2), integrating them using Eq. (fTU|) . and 
finding u by requiring ej(T) = ue/(0). 

The analysis of the eigenvalues of u is simplified by the 
fact that for the equations considered det(u) = 1. This 
can be easily shown by relating u to the 3x3 operator 
U representing the integration of eq. (|10[) in the inertial 
frame. Since Tr(M) = (can be seen directly form the 
expression of the coefficients of M following (fTU)) ). we 
have det(U) = 1 which leads to det(u) = 1. 

The eigenvalue equation for u is x 2 — Tr(u)x + 1 = 
and the problem reduces to finding the value of Tr(u). 
For |Tr(u)| > 2, there is a real eigenvalue x > 1 and the 
orbit is unstable. If |Tr(u)| < 2, both eigenvalues have 
magnitude \x\ = 1 and the orbit is stable. 

For axisymmetric potentials, one cycle takes T = 2tt/(1 
implying that 

Tr(u) = e x + T + e x - T = 2 cosh (2TrVA/n\ , (15) 

where A is given by Eq. ([13]). The condition \Tr(u) > 2| 
is satisfied if and only if A > 0, in agreement with the 
analysis above. 

The resulting values of Tr(u) as a function of inclina- 
tion for the cases of Kozai and oblateness are shown in 
Fig. [TJ As can be seen in the figure, while the Kozai 
case is unstable only above the critical angle, and the 
oblateness case is always stable, both potentials exhibit 
a rich structure with many points satisfying |Tr(u)| = 2. 
For a general axisymmetric potential, |Tr(u)| = 2 occurs 
whenever A < and V— A/f2 is an integer (see Eq. (1151) ). 
As shown below, these points can become unstable once 
small non-axisymmetric potentials are added. 

Consider next the combined Kozai and oblateness po- 
tential. One example is a satellite of the Earth per- 
turbed by its oblateness and the Moon. For an Earth 
orbiting satellite perturbed by the Earth's oblateness 
and the Moon, e is related to the semi-major axis by 
a 6.26i?Earth£ _1 ^ 5 with the secular time scale given 
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FIG. 1: The value of Tr(u) as a function of inclination, given 
by Eq. (|15[1 for the axisymmetric Kozai and Oblateness po- 
tentials. Instability corresponds to |Tr(u) | > 2. The red cross 
corresponds to the Kozai critical inclination. 



-1 





t 




r i 
\ | 






I | 

1 , 

I 

1 1 

1 1 
. 1 1 
1 ' 






I 

1 1 
I 1 
ll 






1 





10 20 30 40 50 60 

i (deg) 

FIG. 2: The value of Tr(u) as a function of inclination for 
combined Kozai and oblateness potential JTJl with e = 0.25 
for iKoz-Obi = (dashed black line) and iKoz-Obi = 23° (red 
solid line). Arrows indicate inclinations used in Fig. [3] 

by i S cc ~ 77e 3 / 10 yr. The value of Tr(u) for the case of 
e = 0.25 is presented in Fig. [21 for ?Koz-Obi = (dashed 
black line) and iKoz-Obi = 23° (red solid line, 23° is the 
angle between Earth's equatorial plane and the ecliptic 
plane.). As can be seen the aligned case has one unstable 
zone above the critical inclination corresponding to Eq. 
([14]) (for e < 1) while the misaligned case has additional 
unstable zones roughly corresponding to the inclinations 
in which |Tr(u)| = 2 in the axisymmetric case. 

The results of numerical integrations of the full (non 
linear) equations of motion Eqs. ([3]) for a few initial incli- 
nations indicated by arrows in Fig. [5] (with initial values 
e = 10~ 4 ,£! = —tt/2,uj = 0.5) are presented in Fig. [3] 
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|Tr(u)| — 2 > accurately captures the regions where e 
reaches significant values. 
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FIG. 3: Evolution of the eccentricity for the combined Kozai 
(dashed red line) and Oblateness (solid black line) potential 
with iKoz-Obi = 23°, for a set of initial inclinations plotted in 
different colors, as indicated by arrows in Fig. [2] 
Unstable (stable) orbits are shown in solid (dashed) lines. 




FIG. 4: The maximum eccentricity as a function of initial 
inclination for the combined Kozai and Oblateness potential 
with JKoz-Obi = 23° (red solid). The value of |7Y(u)| — 2 from 
linear analysis (dashed-dotted green) and the maximum ec- 
centricity obtained in the same time for the aligned combined 
potential (iKoz-Obi = 0, blue solid) are shown for comparison. 



Integrations with initial inclinations in unstable (stable) 
zones are shown in solid (dashed) lines. As can be seen, 
in the unstable zones the eccentricity grows exponentially 
to significant values. 

The maximum eccentricity reached in the interval 
t = — 20i scc as a function of initial inclination (with 
the same initial values e = 10~ 4 , £1 = — 7t/2,cj = 0.5) is 
shown in Fig. [4] (solid red). For comparison, the value of 
|TV(u)| — 2 is shown on the same figure (dashed-dotted 
green). As can be seen, the linear instability condition 
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FIG. 5: Results for combined potential with e = 2.5. The 
upper panel shows the value Tr(u) for aligned (black dashed) 
and iKoz-Obi = 23° (red solid). The bottom panel shows the 
corresponding maximum eccentricity as a function of inclina- 
tion. |Tr(u)| — 2 for JKoz-Obi = 23° is shown for comparison 
(green dashed). 

Similar results for the combined potential with the 
same alignments iKoz-Obi = 23° and e = 2.5 are pre- 
sented in Fig. [5l Similarly, a rich structure of stable and 
unstable zones is obtained. 

Finally, we perform stability analysis on the quadratic 
potential given in Eq.® with a = and 6 = 0.1, which 
corresponds to a small non-axissymmetric component 
added to the Kozai potential (that is, a = 5 = in 
Eq. ([5])). The results are shown in Fig. [B] . As can be 
seen in the upper panel, due to the non-axissymmetry, 
a new instability zone is formed near io ~ 35 deg, where 
Tr(u) = 2 for Kozai. As illustrated in the red solid 
line in the bottom panel, within 40i scc , substantial ec- 
centricities are obtained in this new zone. Interestingly, 
the attainable maximum eccentricities beyond the Kozai 
critical angle are significantly larger than those in the 
Kozai case (black dashed). In fact, upon closer inspec- 
tion, we find maximum eccentricities reach as high as 
~ 0.999 during short episodes in which the inclinations 
cross 90 deg, similar to the behaviors obtained when the 
octupole contribution is added to the Kozai quadrupole 
potential [HI, HU . Detailed discussions are deferred to a 
future publication. 

Upon completion of this work, we learned that the sta- 
bility analysis on the combined potentials of Kozai and 
oblateness was performed in 
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FIG. 6: Results for the quadratic potential with a = and 
8 = 0.1. Tr(u) as a function of initial inclination is shown in 
red solid line in the upper panel. For comparison, the case 
a = 8 = (Kozai) is plotted in black dashed. The bottom 
panel shows the corresponding maximum eccentricity in the 
interval t — — 40t scc . |TV(u)| — 2 for 5 = 0.1 is shown for 
comparison (green dashed). 
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